# Load packages
librarian::shelf(dplyr, estimatr, haven, psych, skimr, tidyr, quiet = TRUE)

# Load the Stata data file
data <- read_dta("Data/Meta/broockmanEtAl_2023_study5.dta")

# Drop participants who didn't finish the survey
data <- data %>% subset(finished_survey == 1)

# Drop participants in the negative_trust_game and outparty_friend conditions
data <- data %>%
    subset(treatment %in% c("control", "lees_cikara", "ahler_sood"))

# Recode variables
data$treatment <- data$treatment %>%
    factor() %>%
    relevel(ref = "control")

data$normswhenpossiblefieldinpartyicp <- 8 - data$normswhenpossiblefieldinpartyicp

# Rename norms items
data <- data %>%
    rename(
        norms_cooperate = normswhenpossiblefieldinpartyicp,
        norms_governor = normsifafieldinpartyicgovernorof,
        norms_court = normsfieldinpartyicelectedoffici,
        norms_journalist = normsifajournalistaccusesafieldi,
        norms_pollplaces = normsfieldinpartysshouldreduceth,
        norms_mail_in_voting = normsmailinvotingshouldbebannedi
    )

# Check if norms load together
data %>%
    select(starts_with("norms"), -norms_scale) %>%
    alpha()

# Generate variables
data$norms_scale <- data %>%
    select(starts_with("norms"), -norms_scale) %>%
    rowMeans(na.rm = TRUE)

data <- data %>%
    mutate(
        affpol = therm_inparty - therm_outparty,
        eliteaffpol = therm_inpartyelites - therm_outpartyelites,
        pid_strength = abs(pid7 - 4),
        pid_republican = ifelse(pid7 %in% 1:3, 0, ifelse(pid7 %in% 5:7, 1, NA)),
        pid_string = ifelse(pid_republican == 0, "D", ifelse(pid_republican == 1, "R", NA))
    )

# Get ranges
max(data$affpol, na.rm = TRUE) - min(data$affpol, na.rm = TRUE)
max(data$eliteaffpol, na.rm = TRUE) - min(data$eliteaffpol, na.rm = TRUE)

max(data$sd_comfortneighbor, na.rm = TRUE) - min(data$sd_comfortneighbor, na.rm = TRUE)
max(data$profile_worker, na.rm = TRUE) - min(data$profile_worker, na.rm = TRUE)
max(data$profile_attractive, na.rm = TRUE) - min(data$profile_attractive, na.rm = TRUE)

max(data$norms_scale, na.rm = TRUE) - min(data$norms_scale, na.rm = TRUE)

# Check correlations
data %>%
    select(
        affpol, eliteaffpol, sd_comfortneighbor, profile_worker,
        profile_attractive, norms_scale
    ) %>%
    cor(use = "pairwise.complete.obs")

# Run regressions
lm_robust(affpol ~ treatment + pid_strength + pid7, data = data)
lm_robust(eliteaffpol ~ treatment + pid_strength + pid7, data = data)

lm_robust(sd_comfortneighbor ~ treatment + pid_strength + pid7, data = data)
lm_robust(profile_worker ~ treatment + pid_strength + pid7, data = data)
lm_robust(profile_attractive ~ treatment + pid_strength + pid7, data = data)

lm_robust(norms_scale ~ treatment + pid_strength + pid7, data = data)
